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Abstract 

We present new integral representations in two dimensions for the elas¬ 
tance problem in electrostatics and the mobility problem in Stokes flow. 
These representations lead to resonance-free Fredholm integral equations 
of the second kind and well conditioned linear systems upon discretiza¬ 
tion. By coupling our integral equations with high order quadrature and 
fast multipole acceleration, large-scale problems can be solved with only 
modest computing resources. We also discuss some applications of these 
boundary value problems in applied physics. 


1 Introduction 

A classical problem in electrostatics is the analysis of capacitance. Briefly stated, 
for an open system in two dimensions, this concerns a collection of N disjoint, 
bounded regions, denoted by Dj, with boundaries T,-, all of which are assumed to 
be perfect conductors. Setting the potential (the voltage) on the jth conductor 
to (j)j for j = 1,... ,N, one would like to determine the net charge qj which 
accumulates on the conductors Dj for j = 1,..., N. Since electrostatics is 
governed by a linear partial differential equation (the Laplace equation), there 
is a matrix, denoted by C, such that 

q = C 0, (1) 

where 4> = ..., 4>n) and q = (qi ,..., <pv)- The matrix C is referred to as 

the capacitance matrix. 

Less well-studied is the converse problem, where a fixed amount of charge 
is placed on each of the N conductors, and the goal is to determine the corre¬ 
sponding, unknown, potentials. The matrix corresponding to that mapping is 
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called the elastance matrix l], denote by P, with 

0 = Pq. 

We should note that the goal is achievable, since the voltage on a perfect con¬ 
ductor is constant from Maxwell’s equations [2], P is the inverse of C, in a 
suitably-defined space, discussed in greater detail in the next section. Given 
either the capacitance or elastance matrix, it is straightforward to compute the 
electrostatic energy E of the system [2j, since 

E= \ (t,T(l= l ( t ,Tc <l>= b qTpq ' 

In some contexts, particularly in chip design, the capacitance problem is more 
typical j3,[4}[5j[6]. In others, including some quantum mechanical settings, the 
elastance problem arises more naturally [ 7,8] 9 . The elastance matrix is some¬ 
times referred to as the charging energy matrix. 

A similar duality exists in problems of Stokes flow. Given N disjoint, rigid 
bodies, denoted by £b, with boundaries Tj, with prescribed translational and 
rotational velocities, (vi,Wi), the resistance problem consists of determing the 
corresponding forces and torques (Fj,T,) on each of the bodies. The mobility 
problem is the reverse; given prescribed forces and torques on each of the rigid 


bodies, find the corresponding velocities (see, for example, 10 11 12 13]). The 
mappings R and M such that 

F = RU and U = MF 

are known as the resistance and mobility tensors. Here, U = (vi ,uii, ■ ■ ■, vn, w/v) 
and F = (FijTi,... ,F N ,T N ). 

Reformulating the governing partial differential equation as a boundary in¬ 
tegral equation is a natural approach for the above problems, since this reduces 
the dimensionality of the problem (discretizing the boundaries T^ alone) and 
permits high order accuracy to be achieved in complicated geometries. More¬ 
over, boundary integral equations can be solved in optimal or nearly optimal 

14,1 5p6p7 


time using suitable fast algorithms 


and satisfy the far field bound¬ 


ary conditions necessary to model an open system without the need for artificial 
truncation of the computational domain. 

In this paper, we will restrict our attention largely to the formulation of 
suitable integral equations for the elastance and mobility problems. There is a 
substantial literature on the development of well-conditioned second-kind Fred¬ 
holm equations to address the capacitance problem, which we do not seek to 
review here. We simply note that to apply C to a vector is equivalent (in the 
two-dimensional setting) to solving the Dirichlet problem: 


Au(x) = 0 xel 2 \ (u^LjA) 


u\ Tj = 


'3 > 


( 2 ) 

(3) 
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together with the radiation condition that u(x) be bounded as |x| —> oo. The 
boundedness of u(x) enforces charge neutrality on the collection of conductors. 
To see this, note that standard multipole estimates imply that 

U (x)^C+£log|x| 

as |x| —> oo, where C is constant and Q is the net charge induced on all the 
conductors. Thus, there would be logarithmic growth of the potential at infinity 
if the system were not charge neutral. 

Remark 1. It is worth noting that the constant C cannot be specified indepen¬ 
dently. If, for example, (j>j were set to 1 on all conductors, the solution to the 
Dirichlet problem must be zt(x) = 1 (under the assumption of charge neutrality), 
so that C = 1. 


The Dirichlet problem, as noted above, is well-known to have a unique solu¬ 
tion and a variety of well-conditioned integral equations have been derived for its 
solution (see [5,18,19,[20j[2l] and the references therein). The resistance prob¬ 
lem involves solving the Stokes equations with boundary conditions imposed 
on the velocity, for which there are, again, a large number of well-conditioned 
formulations 


12 13 22 23 24 


Suitable integral representations have been developed for both the elastance 

or 


and mobility problems, often in the form of first kind integral equations 25 


second kind integral equations with N additional unknowns and N additional 
constraints 26 . While these have been shown to be very effective, when N is 


large and the geometry is complex, it is advantageous to work with formulations 
that are both well-conditioned (formulated as second kind boundary integral 
equations) and free of additional unknowns and constraints. We develop such 
an approach for the electrostatic problem first, in section [2j followed by the 
Stokes mobility problem in section [3] We illustrate their effectiveness with 
numerical examples in section [4] and discuss generalizations in section [5] 


Remark 2. We note that the elastance problem can be interpreted as a special 
case of a modified Dirichlet problem, and second kind Fredholm integral equa¬ 
tions for the modified Dirichlet problem are developed and discussed in 1 20 [ /. In 
our representation, we compute the physical charge density on each conductor 
directly and stably. For the representation discussed in \20f , the charge density 
is given in terms of a hypersingular integral. 

Second kind integral equations for mobility problems without additional con¬ 
straints were developed earlier by Kim and Karrila \2fj , using the Lorentz re¬ 
ciprocal identity. Our derivation, which leads to essentially the same integral 
equation, is direct - based on the physical principle that the interior of a rigid 
body must be stress free. Finally, the mobility problem can also be solved using 
a double layer representation without additional unknowns. For a detailed dis¬ 
cussion of this representation, we refer the reader to Our formulation 

has the advantage that certain derivative quantities, such as fluid stresses, can 
be computed with weakly singular instead of hypersingular kernels. 
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2 The Elastance Problem 


Given the set of conductors {Dj} i=1 , let us assume that the boundaries T, : = dDi 
are positively oriented. We will denote by T the total boundary T = U^ 1 P i 
and by n x the outward normal at x £ T. We let E = R 2 \ (u f =1 Di) denote the 
exterior domain. For the sake of simplicty, we assume that the conductors have 
smooth boundaries (Fig. [l]) 



Figure 1: Three, smooth bounded conductors in the plane, with the exterior 
domain denoted by E. For x on the boundary, n x represents the outward 
normal. 

Application of the elastance matrix to a vector of charge strengths q = 
(qi, q 2 ■ ■ ■ qN) is equivalent to the solution of the following boundary value prob¬ 
lem for the potential u(x) in the exterior domain E: 


Au(x) =0 x e E 

(4) 

e 

n 

II 

(5) 

du 


Tn ds, = qj 

(6) 

u (x) —>• 0 as x — ► oo . 

(7) 


Here, it(x) and the constants {4>j}^ =1 are unknown. As noted in the introduc¬ 
tion, it is a consequence of the Maxwell equations that the potential on each 
distinct conducting surface is constant, so that the boundary condition © cor- 
responds to the physical problem of interest. Q enforces the desired charging 
of the individual conductors, and ([7]) corresponds to setting the potential at 
infinity to zero (ground). 

Remark 3. ('Charge neutrality/- It is often said that the elastance matrix P 
is the inverse of the capacitance matrix C. Unfortunately, it is straightforward to 
verify that, for the vector of potential values <p Q = (1,..., 1), we have C </> 0 = 0, 
so that C is not actually invertible. Likewise, the elastance boundary value 
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problem, as stated above, cannot be solved unless q = (q\,q 2 • • • <7 at) satisfies 

N 

E^' = °- (8) 

3=1 

Otherwise, u(x) would have logarithmic growth at infinity. That is the sense 
in which P is the inverse of C - as a map defined on the space of mean zero 
vectors in M. N . 


To prove uniqueness for the elastance problem, we will need the following 
lemma 


29 


Lemma 1. Let u be a harmonic function in the exterior domain E defined 
above, satisfying the condition ©. Let B r (0) be the ball of radius R centered 
at the origin and let DBr (0) be its boundary. Then, there exist M, Rq such that 
su Ppb h (o) |Vu| < %: for all R > R 0 . 

Lemma 2. (Uniqueness). Suppose that u satisfies equations ©■ © and 
©, with qi = 0 for i = 1,..., N. Then, u(x) = 0 in the exterior domain E. 


Proof. For sufficiently large R, we may write 


0 = f uAudV = f d.Sx — f \SJu\ 2 dV. 

Jedb r(0) JdB R (0) On i=1 Jri on J_Ens H (o) 


Since u(x) takes on some constant value <fi on Fj, we may write 


f |Vu| 2 dV = [ 
JedB r (0 ) J d 


r) N 

OB fl (o) dn 


>dB R ( 0) 


du , 

u—ds x , 
on 


du 

dn 


ds x 


since the qt are all zero. From Lemma[l] the boundedness of u, and the monotone 
convergence theorem, it is easy to see that 


\Vu\~dV 


IE 


lim / |Vu| 2 dF 

R ^°°JEnB R (0) 


lim 

R—> oo 


dB R ( 0) 


du, „ 

u—ds x = 0. 
dn 


Thus, Vu = 0 in E and u must be a constant. From the decay condition at 
infinity, u = 0 as desired. □ 


To develop an integral equation for the elastance problem, we will use the 
language of scattering theory. That is, we will construct an “incident field” 
which satisfies the charging conditions © but not the boundary conditions ©. 
We will then solve for a “scattered” field which forces the conductors to be 
equipotential surfaces without changing the net charge on any of the r^. (This 
will also yield a proof of existence of solutions.) 
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Remark 4. In physical terms, one can think of the problem as follows: imagine 
that we simply deposit charge uniformly on each conducting surface to satisfy 
the charging condition. This will be our incident field. The charges will then 
redistribute themselves on each T* so that they are equipotential surfaces. The 
total field will be defined by that new, equilibrated charge distribution. 


2.1 Mathematical preliminaries 

Let 7 be a smooth closed curve in R 2 and let D T denote the domains corre¬ 
sponding to the interior and exterior of 7 . Let n x be the unit outward normal 
to the curve 7 and let n 0 = n Xo for x 0 £ 7 . Let p : 7 — > R be a continuous 
function. The single layer potential is defined by 

SjH (x) = f G(x,y)p(y)ds y , (9) 

J1 

where G(x, y) is the fundamental solution for the Laplace equation in free space: 


G(x,y) = - — log|x-y| . 


( 10 ) 


Lemma 3. |20j29|30j Let Sj/Jj (x) be a single layer potential with charge density 
p defined on 7 . Then S 1 p (x) is harmonic in R 2 \7 and continuous in R 2 . The 
single layer potential satisfies the jump relations 


lim 

x->x 0 

xea* 


dS 1 p (x) 
<9n 0 ,± 


~F 2 A 1 ( x o) 


dG (x 0 , y) 
d n 0 


p (y) ds y 


( 11 ) 


where j>' indicates the principal value integral over the curve 7 and the subscripts 
— and + denote the limits of the integral from the interior and exterior side, 
respectively. Furthermore, 


[ <9S> (x) 


p (x) ds x 


l dS 1 p (x) 


For a closed curve u j C D + , we also have 


dS 7 p (x) 

drijr 


ds x = 0. 


( 12 ) 


(13) 


Finally, 


where 


S'yP ( X ) + ^-Qlog (x) 


0 as |x| —> 00 , 


Q= T(y)ds y . 

J -y 


(14) 
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The double layer potential D 7 p{x) is the potential due to a surface density 
of dipole sources on 7 , aligned in the normal direction to the curve: 


f dG{x, y) 

D 7 /i (x) = j ——- \x (y) ds y 


(15) 


Lemma 4. 


is harmonic in 


I 


,25 3Q| / Let Z) 7 /r (x) 6 e a double layer potential. Then, D^p (x) 
C 7 and satisfies the jump relations: 


x hm f) £> 7 M = ±~p (x 0 ) + 

X6D* 


dG (x 0 , y) 
dn y 


p (y) ds y . 


Furthermore, 


f 3G(x,y) 

7 

dG(x,y) 

dn v 


ds x 


ds x 


l-i xef)- 

|o xgd+ 

XG7 ’ 


(16) 


(17) 

(18) 


and 

\D 1 p (x)| —>• 0 as |x| —> 00 . (19) 


2.2 Charging the boundaries with an incident field 

In the elastance problem, perhaps the simplest way to allocate the net charge 
qt to the boundary T* is to define a constant charge density 

0i( x ) = ]|rj-, (20) 

for x on the curve C, where |Tj| denotes its length. We can then define cr (x) = 
(cti (x) ,a 2 (x)... a N (x)) and 


Uinc (x) = Srcr (x) 
where Sr is the operator given by 

n N r 

S r cr (x) = s r : a 3 ( x ) = / G ( x ’ y ) ( y ) ds y ' 

j =1 J = 1 


( 21 ) 


From (12) and (13), we have 


L dUl dn ] ds x - L ^ (x) dSx “ 1C I L dSx - q ° 


( 22 ) 


for j = 1, 2,..., TV. Thus, Ui nc satisfies the charge constraints 
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2.3 The scattered field 

We now seek a scattered field 


(23) 


u ac (x) = Sr A* (x) 

such that u(x) = u inc (x) + it sc (x), where /r(x) = (^(x), /U 2 (x),..., /zjv(x)) 
with /Tj an unknown charge density on the boundaries 1^-. To ensure that no 
additional net charge has been introduced on any of the conductors, we impose 
the N integral constraints on /x (x): 

/ fjLj (x) ds x = 0 . 

•'U 

If we can find such functions /ij(x), then 

u(x) = u inc (x) + u sc (x) = S r (n + cr)(x), (24) 

solves the elastance problem. Physically, (fj,j + cjj ) (x) is the final charge density 
on Fj once the total charge place on the boundary has equilibrated to enforce 
the perfect conductor boundary condition <§• 


2.4 Formulation as a Neumann problem 

Letting u(x) = Ui nc (x) + u sc (x), note first that the corresponding potential is 
also defined inside each conductor. Rather than imposing the boundary condi¬ 
tion ([5]) from the exterior, however, we can make use of the fact that the electric 
field inside each conductor given by V«(x) must be identically zero. Thus, we 
may impose the interior Neumann boundary conditions 


du 
dn- 


(x) = 0 


(25) 


for x G T. Using (11), we obtain the following second kind integral equation: 

T 


I + K ) n (x) = - ( \l + K ) a (x) 


(26) 


for x G T, subject to the constraints 



AU ( x ) ds x = 0, 


(27) 




for j = 1,2,... IV. Here, 


( h 


I = 


h 


\ 


V 

( 

K = 


K 1,1 
K 2 ,i 


In ) 

Kx t2 ... K\ n \ 

A'2,2 ■ • ■ K 2 ,N 


(28) 


(29) 


\ HNA Kn,2 Kn,N ) 

li : C°’ a (Ti) —> C°’ a (r,) is the identity map, Kij : C°’ a (1^) — > C°’ a (T*), and 
Kii : C°’ a (T*) — > C 0,a (fj) are the operators given by 


II 

b 

'^3 

£ 

r dG (x,y) ^ 

Jt 3 9n x 

(y) ds y 

x e r, 

(30) 

Ki.iV = 

f SG^y)^ 

JTi 

(y) ds y 

xeFj, 

(31) 


where C 0,a (r) is the Holder space with exponent a and a > 0. For a related 
treatment of the capacitance problem, see [5]. 


Theorem 1. Letu(x.) be defined as in (J24J). If p(x) solves equations (26) and 
(27), then u (x) solves the elastance problem. 


Proof. Note first that u (x) is harmonic in E by construction. Using (12) and 
(13), the choice of a in equation ([20| and the constraints (27), we see that u (x) 
satisfies the charge constraints ([6j). Furthermore, from equations (| 8 ]| and ( fl4| ), 
it follows that u (x) —> 0 as |x| —> 00 . Since u (x) is harmonic in Di and satisfies 

du 
dn- 


= 0 on Ti, u (x) = Ci for some constant Cj inside Di. By the continuity of 


the single layer potential, u = Ci from the exterior side of Tj as well. 


□ 


Remark 5. (The adjoint operator/- The operator K in equation (26) is a 


compact operator for smooth T. Hence, 1/ + I\ is a Fredholm integral equation 
of the second kind. To study existence of solutions to (|/ + AT) p = /, we shall 
study existence of solutions for the adjoint problem [II + K *) p = f instead, 
where 


K* 


/ K h 

K* 

^ 1,2 

K* 

* 2,1 

K* 

^ 2,2 

K* 
1X 2 ,N 

\ *jV ,1 

K* 

^N ,2 

K * 

1 X N,N 


\ 

/ 


(32) 
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Here, K*^ : C°’ a (T,-) -► C°’ a (r*) and KT : C°-“ (r 4 ) -»■ C°-“ (r 4 ) are the 
operators given by 


K tj a = 


KUa = 


f 0 G(x, y) 
r . 5n y 
9G (x, y) 
dn v 


a (y) dsj 


er (y) dsj 


(33) 

(34) 


It is straightforward to verify that K*jCr(x) = D r .a(x) for x £ Fj. 


Let a = {cFi (x )}^ 1 where each tr* (x), supported on T.;, 
using ,(|18[) and ( |32| ), we may conclude that (,(I + K* 


is constant. Then, 
i j .l ) a = 0. Thus, the 
dimension of the null space of 4/ + K* is at least N. In fact, it is well-known 
that the dimension of the null space is exactly N \20\ 29j. 


Remark 6 . 1/ + K* is the integral operator one would obtain in seeking to 
impose Dirichlet boundary conditions with the potential represented as a double 
layer potential. The double layer potential operator for the exterior, however, is 
range deficient. It cannot represent a harmonic function u (x) in the exterior 
which is generated by net charge in any of the domains D j. To see this, note 
that the net charge is — f r ^ from ( | 12 [ ) , but that the double layer potentials 
satisfies f r = 0 for i = 1,2,... N from ( fl7| ). 


2.5 Existence of solutions 

From the preceding discussion (the existence of a nontrivial nulllspace), it follows 
from the Fredholm alternative that (4/ + /F)/i = / has an N dimensional space 
of solutions, so long as / is in the range of the operator \l + I\. Using our 


representation for the elastance problem, the right hand side in equation (261 
is certainly in the range of the operator 4/ 4 - K. The role of the additional 
N integral constraints is, therefore, to pick out the unique one which doesn’t 
alter the net charge on the N conductors. However, we do not wish to solve an 
overdetermined (non-square) linear system. If we simply discretize the integral 
equation using, say a Nystrom method, with M points on T, we would have to 
solve an (M + N) x M linear system to obtain the desired solution. Instead, we 
propose to solve the integral equation 


1 


/b(x) 


N 

E 

3 =1 


fAj (> 


N 


Mi(x)ds x = (x) - y^FQjo-j(x) 

3 =1 


for x £ Ti, or 


-I 

2 


K + L) n = - 


-I + K } a 


(35) 
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where 


L = 


(U 


\ 


(36) 


Ln ) 


with Li : C°’ a (T,) —> C°’ a (Tj) defined by -Lj/Zj(x) = f r /j, i (y)ds y . 


The following lemma shows that solving (35) is equivalent to solving (26) 


with constraints (27). 


Lemma 5. If h solves equation (35), then /i solves equations (26) and (27) 


Proof. Using equation (13), we observe that f r Ki j/ij (x) ds x = 0 for j i. 
Furthermore, switching the order of integration in Ki^Hi (x) ds x and using 


property (181 of the double layer potential, we see that 


<j> Ki^Hi (x) ds x = 


dG (x, y) 


r, Jr, 


dn x 


Hi (y) dsydsy 


= J Vi(y)£ 


dG (x,y) 
dn x 


ds x ds y 


lH{y)dSy. 


Integrating expression (35) on Tj, we may conclude that 

|Fi| J Hi (x)dsx = 0. 


(37) 

(38) 

(39) 

(40) 


Thus, LiHi{x) = 0, which implies that h satisfies the integral constraints (27) 
and that 


2 1 + K + l) H — (7^1 + k \ h — — (^1 + K ) o' ■ 


(41) 


□ 


Remark 7. For further discussion of the solution of consistent linear systems 
with constraints in the finite dimensional case, see [31 /. 


The following lemma shows that the operator ^1 + K + L has no null space. 
Lemma 6. The operator + K + L is injective. 

Proof. Let h € AT (|l + K + L), i.e. it solves [\l + K + L) h = 0. Following 
the proof of Lemmata we conclude that L/z = 0 and therefore (^1 + K) h = 0. 
Let u = SrH- From the properties of the single layer potential 




(42) 
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By uniqueness of solutions to interior Neumann problem, we conclude that u is 
a constant on each boundary component. Thus, u solves the Elastance problem 
with q.i = 0, as Lp = 0. By uniqueness of solutions to the Elastance problem, 
we conclude that it = 0 in E. Hence, = 0. From the properties of the 
single layer, 



du 
dn + 


= 0 . 


(43) 


Therefore, J\f [\l + K + L) = {0}. 


□ 


By the Fredholm alternative, we conclude that (35) has a unique solution pt. 


3 The mobility problem 

Supose now that we have N rigid bodies immersed in an incompressible Stoke¬ 
sian fluid in R 2 . Let F, , T.; denote the force and torque exerted on rigid body 
Di in a fluid which is otherwise assumed to be at rest and let Vi,LUi be the 
corresponding rigid body motion, where otj is the angular velocity about the 
centroid of Di. The mobility matrix M £ ^ 3 ^x 3 ^ j g jq le p near mapping from 
the forces and torques on the rigid bodies to the respective rigid body motions: 


U = MF 

where U = (vi,wi,. . .,v N ,u N ) and F = (Fi, Ti,... ,F N ,T N ). 

Referring to Fig. 0, let Di now represent the rigid bodies and let E repre¬ 
sent the Stokesian fluid with viscosity /i = \. Further, let us assume that there 
are no other volume forces on the fluid. Let u(x) = (iti(x), ii 2 (x)) represent the 
fluid velocity in E and let (Fi,Ti,.. ., F n,Tn) be the force and torque exerted 
on the rigid bodies. Let x\ = pb f r xds x be the centroid of T*. Let p be the 
fluid pressure and let er be the stress tensor associated with the flow: 

= -P 5 a + (|^T + 55^ ) = -rfa + 2e ( u ) ( 44 ) 

where d l3 is the Kronecker delta, 

e (u) = - [Du + Du r ) (45) 

is the strain tensor associated with the flow, and Du is the gradient of u. 

On the surface of rigid bodies T i; 


f = <r ■ n = 

cu 

Cl2 


ni 


&21 

0"22 


n 2 


(46) 


represents the surface force or surface traction exerted by the fluid on Di , where 


n is the outward normal to Lj. For notational convenience, let x^ = 


-X-2 

x\ 


12 









and V 1 - = 

r _jl i 

dx2 

d 

. Then u (x) solves 

12 

13 


dxi 





—Au + Vp = 0 in E 

(47) 

V • u = 0 in E 

(48) 

u ( x ) lr» = V* + Wi (x - X?) x 

(49) 

j f ds x = J a.nds x = ~Fi 

(50) 

J (f, ( x - xj)- 1 -) ds x = —Ti 

(51) 

u (x) —»• 0 as x —> oo , 

(52) 


where (a, b) represents the Euclidean inner product for vectors a, b G R 2 . Equa¬ 
tions (47) and (48) are the governing equations for Stokes flow in domain E. 
Equation (49) enforces a rigid body motion on where v, ; , w,; are unknown. 
Equations (50) and (51) state the net applied forces and torques are given by 
the known quantities (F i,Ti). Finally, (52) states that the fluid is at rest in the 
absense of forcing. As a consequence of Stokes paradox, there might not exist a 
solution to the set of equations described above. In fact, it can be shown that 


u (x) = O 


N 


- F* log |x| 


Thus, a necessary condition for a solution to exist is that 


(53) 


E F * = °- 


2=1 


(54) 


From 12 13 , it turns out that (541 is also sufficient for a solution satisyfing 


equation (52). To prove uniqueness for the mobility problem, we need the 


following lemmas which can be found in 32 


Lemma 7. If h is a bounded harmonic function in E and n is an integer greater 
than 0, with h = O ( r~ n ) as r oo, then 


OO 

h{r,8) = Y r ~ ka k ( 0 ) ) (55) 

k—n 


which converges uniformly outside Br (0) for some R. 

Let u> be the vorticity corresponding to the flow, defined by 

w=(V- L ,u). (56) 

Lemma 8. If u satisfies equation (1521), then u> = O (r _1 ) asr-> oo. 
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Lemma 9. If uj = O (r ") as r —> oo for integer n > 0, then p = O (r ") as 
r —>• oo. 


Using these two lemmas, it follows that 
Lemma 10. If u satisfies equation ( |52| ), then on dB r (0), 

p (u, n) — to (u 1 , n) = O (r^ 1 ) as r —> oo . 


(57) 


Lemma 11. If u satisfies equations ( |49| ), (|50j) and (51) withFi = 0 andTi = 0, 
then 

J (u, f) ds x = 0. (58) 

Proof. 

J (u, f) ds x = J (vi (x-xj)- 1 ,f) ds x = - (v^Fj) - WjTj = 0. (59) 

□ 


Lemma 12. If u satisfies equations (49), (50) and (51) withFi = 0 andTi = 0, 
then 

(60) 


J p(u, n) - w (u x ,n) ds x = -4a ;, 2 |A| 


Proof. On Ti, e (u) = 0 and ui = 2cc,;. Using Lemma 11 and the divergence 
theorem 


- J (u, f) + 2 u>i (v x + Ui (x - x?) , n) ds x = -4 uf |A| 


(61) 


□ 


Lemma 13 (adapted from 32). If u (x) satisifi.es equations (47), (48), (49), 
(50), (51) and (52) with Fj = 0 and Ti = 0, then 


lim 

R—> oo 


[ (u, f) ds x —?> 0 . 

JdB R ( 0) 


(62) 


Proof. For large enough R , Lemma [l2| yields 
„ N 


UJ‘ 


/ EDBr(0 ) 


'dV = ^ / p (u, n) — w (u x , n) 


i =1 ' 


~f P ( u , n) — w (u x , n) ds x (63) 

JdB R ( 0) 

N 

=-4VW, 2 IAI - / P (u, n) — w (u x , n) ds x . (64) 

JdB R (o) 
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Using Lemma 10 we conclude that 


/ u 2 d.V < oo . (65) 

J E 

Using Lemma [7J we know that 

u(r, 6) = r~ 1 a\{9) + O (r~ 2 ) . (66) 

Integrating w 2 in the annulus B = B r (0) D Bp (0) C , we get 

J uj 2 dV = log(^) J a\(0)d6 + O (r -2 ) . (67) 

Since f B cu 2 dV is bounded, we conclude that a\ = 0 and that w = O (?'“ 2 ). 
Using Lemma [9j we conclude that p = O (r~ 2 ). Thus 


/ [— p(u, n) +w (u x ,n)] ds x —> 0 as R —> 00 . 

JdB R (0) 


( 68 ) 


From equation (64), it follows that w = 0 in E. Using equation (48), we conclude 
that Au = 0 in E. Using the estimate for p and Lemma 111 we get f = O (j"~ 2 ). 
Using this estimate and the decay condition in u at 00 , the result follows. □ 

The following lemma is a modification of the standard proof of uniqueness 


for Stokes flow 12 13 


Lemma 14. If u (x) satisifies equations (47), (48), (49), (50l, (51) and (52) 
with Fj = 0 and Tj = 0 , then u (x) = 0 . 

Proof. Let (•, •) : R 2x2 x R 2x2 be the Frobenius inner product. For large enough 

f?, 


f (e (u), e (u)) dV = f (Du, e (u)) dV 

JEnBji(O) JEC\Br(0) 

(u, e (u) ■ n) ds x — - / (u,Au) 

^ JEnBft(O) 

(u,e(u) • n)ds x - - f (u ,Vp) 

" JEnB r ( 0 ) 


EnB r(0) 

d(EnB R (0 )) 

f 

d(EDB R (0)) 


dV 

dV 


9(EnB s (0)) 

N 


U, -p 


1 0 
0 1 


= _ 9 + n 


OB r (0) 


EnB R (0 ) 

+ 2 e (u)'] .n ) ds x 
(u, f )ds x 


dB R ( 0) 


(u, f) ds x 
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using (48) and Lemma 11 Taking the limit as R —» oo in the above expression 
and using equation (62), we get 


'(u) = 


0 0 
0 0 


x G E. 


Thus, u is a rigid body motion. However since u(x) —> 0 as |x| 
conclude that u = 0. 


(69) 

>, we 
□ 


We construct an integral representation for the mobility problem by direct 
analogy with the elastance problem, with the velocity u(x) playing the role of 
the potential and surface traction f playing the role of charge in the elastance 
problem. (A rigid body has no interior strain or stress, with all the stress 
residing on the surface.) We first construct an “incident” field which satisfies 
the net force and torque conditions on each rigid body but which does not 
correspond to a rigid body motion. We then find a “scattered” velocity induced 


by an additional force vector /i so that the total velocity will satisfy (49) but 


does not change the net force and torque. As in the elastance problem, this can 
be thought of as a redistribution of the surface force. 


3.1 Mathematical preliminaries 

In the remainder of this paper, we will use the Einstein summation convention. 
As above, we let 7 be a smooth closed curve in R 2 and we let D ^ denote the 
domains corresponding to the interior and exterior of 7 . n x = (n x , 1 , n x , 2 ) will be 
used to denoted the unit outward normal at x G 7 and no = (no,i, no, 2 ) to denote 
the unit outward normal at x 0 G 7 . We let £t( x ) = (/xi(x),/z 2 (x)) : 7 —> R 2 be 
a continuous function. 


Following the treatment of 12 13 , the fundamental solution to the Stokes 


equations (the Stokeslet) in free space is given by 


Gm(x, y) = — 


log |x — y| 5tj + 


(xj - yj) (xj - yj) 

|x-y| 2 


MG 1,2. (70) 


The Stokeslet allows us to express the velocity field u = (■u 1 ,u 2 ) induced by a 
point force f = (/ 1 , / 2 ) in the form 


Ui = Gij (x, y) fj . 


(71) 


The single layer Stokes potential is the velocity induced by a surface force on a 
boundary 7 : 



Gi,j (x, y) nj (y) ds y 


for i = 1 , 2 . 


(72) 


Lemma 15. Tef 6> 7 /i(x) denote a single layer Stokes potential of the form (72). 
Then, S 1 n (x) satisfies the Stokes equations in R 2 \y and S-y[i(x.) is continuous 
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in M 2 . Moreover, if we let f(xo) denote the surface traction on 7 corresponding 
to the velocity field S 1 p (x), then 


lim /j,± (x 0 ) = =F \pi (x 0 ) + n 0 ,fc f T, (x 0 , y) pj (y) ds y 


X^X 0 

x.eD d 


where Tij,k (x, y) is the stresslet corresponding to the flow given by 
m , x _ 1 Oi - Vi) {Xj - Vj) (x fc - y k ) 

J- i,j,k ( x ) yj — | 14 

* x-y 


(73) 


(74) 


The notation <f is used, as above, to denote the principal value integral. The 
net force and torque on the domain are given by 


f + ds x = - / n (x) ds x 


f_ ds x = 0 


(75) 


and 


(x-x c ) ,f ds x = — (x-x c ) ,fi ds x , 


(x-x 0 ) -1 ,^ ds x = 0 . 


If to is a closed curve in D + , then 


f ds x = 0 


(x-x c ) ,f) ds x = 0 . 


Finally, 


as |x| —► 00 , where 


1 


1 0 ' 

R 

r 

5 7 /x(x) + — 

log (x) 

0 1 

M 2 . 

/ V(y)dSy 
J 7 


R = 


Xf X1X2 

X 1 X 2 x\ 


(76) 

(77) 

(78) 

(79) 

(80) 

( 81 ) 


The double layer Stokes potential is the velocity field due to a surface density 
of stresslets on the curve: 


'D ( x )j — / T(y, x) p,j (y) dsy . 

J -V 


(82) 


Lemma 16. Let (x) denote a double layer Stokes potential of the form 
(82). Then, 2? 7 /x(x) satisfies the Stokes equation in R 2 \7 and the jump rela¬ 


tions: 


Hm 'D 1 p i = zfci/ij (x 0 ) + T j4tk (y, x 0 ) pj (y) n y , k ds y . 


X—>-X 0 

xeD ± 


(83) 
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Furthermore, 


TlydSy = 

■os 

1 o 

x e D 
x G E ’ 

(84) 

Thy dSy - 

Sij 

2 

x e 7 . 

(85) 


Letting euj be the standard Levi-Civita symbol, 
j EilmyiTm,j ,k (y- x) kly k dSy = ^ 
^ ^ilmyi^mj.k (y• x) U.yA, dSy — 


Finally, 


I'Djfi (x)| —> 0 as |x| —► oo . 


—cujXi xefl 

0 x G £ ’ 

(86) 

^ilj Xl 

2 

(87) 

1 — )• OO . 

(88) 


3.2 Applying the net force and torque as an incident field 


We construct a velocity field uj nc (x) in the exterior domain E, due to set of 
surface force densities {pj (x )} i ._ 1 on the boundaries {Fj}^ 1 , which satisfies 
the force and torque constraints (50) and (51). Each pj is a vector density 
Pj = Letting p (x) = (p 1A (x), p 2 ,i (x)... pi, N (x), p 2 ,JV (x)), we 

define 


u inc (x) = Srp (x) , 
where <Sr is the operator given by 


(89) 


n N r 

s rP(^)i = ^2^r d P j {x) i = ^2 G i,k(^y)Pkj(y)ds y * = 1,2. (90) 

j =i j =l ,Jr j 

If we now let f ? denote the surface force on Tj corresponding to the velocity 


field u i nc and make use of equations (751 and (78), we obtain 


F i = - / f jds x = Pj (x) ds x for j = 1,2, ...N. (91) 

J r,- Jr, 


Using equations (76), and (79), we obtain 


T 3 = - l (( x - x j)' L : f j) ds * = ((x-X^^^^dSx. (92) 
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Thus, any choice of pj (x) which satisfies equations (91) and (92) will define an 
incident field that enforces the desired force and torque conditions. We will use 
the simple formula 




\ T j\ 

where |Tj| is the length of and Wj = f r 


W, 


(93) 


‘ ds* 


3.3 The scattered field 

We now seek a “scattered” velocity held u sc (x) induced by unknown surface 
force densities { p ;j (x)} _ ± on the boundaries {Tj}” =1 . Each pj is a vector 
density pj = (pij, P 2 ,j)- These densities correspond to a redistribution of 
surface forces that will be used to enforce the rigid body boundary conditions 
without affecting the net force and torque. We let 

p (x) = (p lt 1 (x) , P2,l (x) . . . Pi ,JV (x) , p 2 ,N (x)) 


and dehne 


u sc (x) = S r p (x) . (94) 

To ensure that no additional net forces or torques are introduced on the surfaces 
Tj, we need to impose 3N integral constraints on p (x), namely 



The total velocity held is given by 

u (x) = \i inc (x) + u sc (x) = Sr ( p (x) + p (x)) . (97) 


3.4 Reformulation as an interior boundary value problem 

The function u(x) = 5r(A*(x)+p(x)) also represents the velocity held inside the 
rigid bodies. Since there is no internal stress in a rigid body, the stress tensor 
er must be identically zero within Di. Thus we will seek to impose 


f_ = (<r ■ n)_ = 0 


(98) 


for x £ T. Using equation (73), we obtain the following Fredholm integral 
equation of the second kind: 


^1 + /C ) p (x) = — ( ^1 + 1C ) p (x) x £ T 


(99) 
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which we subject to the constraints 


/ fiij (x) ds x = 0 , 

■' r 3 

j (( X ~ ds x = 0, 


( 100 ) 

( 101 ) 


where 




I = 


V 


and 


JC = 


( £1,1 /Ci, 2 

/C2,l A^2,2 

V /Cjv,1 /Cjv,2 


Ijv / 

/Ci.iv ^ 

&2,N 

/Cat,at 7 


Here, I; : C 0,a (F;) x C 0, “ (r^) — » (7°’“ (r*) x (7°’“ (r*) is the identity map, and 
JCij : c°’ a (Tj) x C°’ Q (r^) ->• c°’ a (r») x c°’ a (Ti) is the operator given by 


for i ^ j and 


(fci,jP)k = n *,i / T k,m,l ( x , y) Pm (y) dsy 

Jr, 

(&i,iP) k = «x,i T fe , m , z (x, y) p m (y) ds y 


x G Tj 


xer, 


Theorem 2. Let u (x) &e t/ie total velocity, defined in (97). If p (x) solves 


equation (99), together with the constraints (100) and (101), then u (x) solves 


the mobility problem. 

Proof, u (x) clearly satisfies the Stokes equations in E by construction. Using 


equations (75), (76), (78) and (79), the choice of p in equation (93), and the 


constraints ( 100 ) and ( 1011 , we see that u (x) satisfies the net force and torque 


conditions (50) and (51). Furthermore, from (54) and (80), it follows that 
|u(x)| -G 0 as |x| —► oo. Since u (x) solves the Stokes equations in D, ; and 
satisfies f_ = 0 on Tj, u must be a rigid body motion. By the continuity of the 
single layer potential, u must define a rigid body motion from the exterior as 
well. □ 


3.5 Existence of solutions 

It is well-known that |l + /C has a 31V-dimensional null space 
from the Fredholm alternative that (^1 + K.) p = g has an 3 N dimensional space 


13 


It follows 


of solutions, so long as g is in the range of the operator \l + K. From (99), this 
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is clearly the case, and the purpose of the 3 N integral constraints is to select 
the particular solution that doesn’t alter the net forces and torques. As for the 
elastance problem, however, we do not wish to solve a rectangular linear system. 
If we discretize the integral equation using Nystrom quadrature, with M points 
on r, we would have to solve a (2 M + 3 N) x 2 M linear system to obtain the 
desired solution. Instead, we propose to solve the integral equation 


( x ) + KijVj ( x ) + / Mi (y) ds 
2 i= i J?i 

+ (x-x c )- L ^ ((y-x^,, 




= ^ \ Pi ( x ) - p o W x G r i 


or 


where 


3 =i 


^I + /C + LU = - Ql + JC)p, 


( 102 ) 


(103) 


(U 


L = 


V 


\ 

L w ) 


(104) 


with L, defined by 

L iMi ( x ) = J r (y) ds y + ( x - x i)^ ((y - *i) ± , 


Mi 


ds y . (105) 


The following lemma shows that solving (1031 is equivalent to solving (991 


with the constraints (1001 and (1011. 


Lemma 17. If fi solves (103), then it solves (99), (100) and (101) . 


Proof. Using equation (861, we see that f r ^(x — x^) -1 ,/Cjj/L ij (x)^ ds x = 0. 
Similarly, using equation (87), we see that 


j r (( x - x i)' L ,^i,iMi( x ))^x = -^ ((x-Xjf^ifx^dSx. (106) 
Since xf is the centroid of T.j, f r (x — x?) -1 ds x = 0. Taking the inner product of 


(103) with (x — x^) -1 , integrating the expression over T^, and using the equations 


above, we obtain 



(107) 
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From (84), we observe that f r /C,;j/x 5 (x) ds x = 0 for j ^ i. Furthermore, 
switching the order of integration in f r (x) ds x and using property (85) 

of the double layer potential, we find that 


J (x) ds x = -^ j n, (x) 


ds x 


(108) 


Integrating the expression (103) on T, and using the fact that 


we may conclude that 


J r ((y - x-)' L , ^ (y)) ds y = 0, 
\ Ti \ f Mi(x)ds x = 0. 


(109) 


Thus, /x satisfies the integral constraints (100) and (101), implying that 
Lj/Xj (x) = 0 and that 


1 . 


1 . 


I + /C + L /x= -I + /CWx = — -I + /C p 


1 . 


( 110 ) 


□ 


The following lemma shows that the operator |l + /C + L has no null space. 
Lemma 18. The operator 11 + K, + L is injective. 

Proof. Let /x £ Af(^I + /C + L), i.e. it solves (|l + JC + L) /x = 0. Following 
the proof of Lemma we conclude that /x satisfies the force and torque con¬ 
straints given by equations ( 100| ) and (101). Thus L/x = 0 and (|I + /C) /x = 0. 
Let u = <Sr/x. Let f_ and f + denote the interior and exterior limits of the surface 
traction corresponding to the velocity field u, respectively. From the properties 
of the Stokes single layer potential 


f_ = ( ^1 + AC ) n = 0. 


(HD 


By uniqueness of solutions to interior surface traction problem, we conclude 
that u is a rigid body motion on each boundary component. Thus, u solves 
the mobility problem with F z = 0 and Xi = 0. By uniqueness of solutions to 
the mobility problem, we conclude that u = 0 in E. Hence, f+ = 0. From the 
properties of the Stokes single layer, 


/x = f _ f + = 0 . 
Therefore, Af Ql + AC + L) = {0}. 


( 112 ) 

□ 


By the Fredholm alternative, therefore, (103) has a unique solution /x. 
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4 Numerical Examples 

The fact that the capacitance and elastance problems are inverses of each other, 
and that completely different techniques can be used for their solution, permits 
a robust test of the performance of our method in arbitrary geometry (without 
an exact reference solution). The same is true for the resistance and mobility 
problems. 


4.1 The elastance problem 

Suppose now that we solve the capacitance problem discussed in section [2] using 
known techniques (see, 20 , for example). That is, given prescribed potentials 


ct>j on a collection of perfect conductors with boundaries Tj, we may obtain the 
charges induced on each conductor. We can then solve the elastance problem 
with these charges as input, using the representation in section[2j and verify that 
the corresponding potentials are those used in the original capacitance problem 
setup. We emphasize that the integral equations used for the capacitance and 
elastance problems are not inverses of each other, so this provides a nontrivial 
test of accuracy. 

More precisely, following the discussion in section [2j we consider the domain 
exterior to N perfect conductors Di whose boundaries are given by T^. We pre¬ 
scribe potentials <f>i on the boundaries T^ and solve the capacitance problem to 
obtain the net charge qi on T^ and also the potential at oo, Uoo = limix^^ u (x). 

We use these charges as input for the elastance problem and to compute 
the potentials induced on the conductors, letting (T l e i denote the uniformly 
distributed charge defined in terms of qi, as in section [2(2} /ii >e q as before, 
represents the unknown density on T^ for the elastance problem and 


w(x) = U inc (x.) + U sc (x) + U ao = SriVel + <7ei)( X ) + u oc ■ 


We then solve 


-I + K + L J Hei — 


2 ^ + K ) 


(113) 


(114) 


where u 00 is the potential at oo computed in the capacitance problem, and the 
operators K and L are described in section |2.5| This is a small modification 
of the representation presented in section |2.2| to account for the potential at 
oo. After solving for we may check the accuracy with which u in equation 
113 equals the potential <p 3 on Tj, the potentials prescribed in the original 


capacitance problem. 


4.1.1 Two disc test 

We first consider the case of two unit discs separated by a distance d. This is 
useful because the exact solution is known and because we wish to study the 
physical ill-conditioning of the problem as d —> 0. 
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Figure 2: Discretization of the discs for Elastance example. 

In the context of Fig. [2j we set u|r 4 = 4>i for i = 1,2 (where i = 1 cor¬ 
responds to the left disc), with <j)\ = 0.209 and 0 2 = —0.123. We consider 
d = 0.5,0.05,0.005. 

We use a Nystrom discretization, based on subdivision of the boundary 
into panels, with Gauss-Legendre nodes given on each panel. Let denote 
the jth node on the ith panel on boundary component l. Let ay j.q e (. i-ii.jj.ei 
denote the density evaluated at We use a recently developed quadrature 
scheme, denoted by GLQBX (global + local quadrature by expansion) [33,34] 
for evaluating the layer potential K in equation (114). This scheme is a robust 


extension of the QBX method of 135], guaranteed to yield high order accuracy 
even when boundaries are close-to-touching. We use an iterative GMRES-based 
solver to obtain to obtain jj, e i, and iterate to a relative residue of 10 -6 . As d —> 0, 
the problem becomes physically ill-conditioned, requiring an increasing number 
of iterations. To improve the rate of convergence, we use an L 2 -based rescaling 
of the unknowns 


36 


. That is, 


we use n s e ^ ale = 

the discrete 2-norm approximates the L 2 norm where r,; is the length of panel i. 


ri as unknowns, so that 


We iterate the following discretized linear system 


D ( X -I + K + L 


D 


-1 scale 
H-el 




K 


Oel 5 


(115) 


where K and L are discretized versions of K and L respectively, D is the diagonal 
operator described given by D^ e i = nlf ale . 

In Fig. [3] we plot the net charge density ai, e i + fii >e i for the three different 
values of d. In Fig. [4] we plot the potential using the off-surface evaluation 
method of [37 , whose development initially led to QBX. (The option of off- 
surface evaluation has been incorporated into our QBX software.) 

From symmetry considerations, 02 , e ; + H 2 ,ei = ~ (ay, e ; + Hi, e i )■ 
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• • d=0.5 

• • d=0.05 

• • d=0.005 . 




11.1. 



4 -3 -2 -1 


2 3 4 


Figure 3: Solution of integral equation for the elastance problem, o \+ ni, e i 
as a function of d. 



Figure 4: Contour plot of u in the exterior of the two discs for fa = 0.209 and 
fa = —0.123 for d= 0.05. 

As noted earlier, the two disc Dirichlet problem has an analytic solution. 
For this, suppose that the left disc is centered at xf = (—1 — f,0), that the 
right disc at x§ = (l + |, 0), and that the discs are held at constant potentials 
fa and fa. Then, the exterior potential is given by 

“.,(x)--giog( !^j°;°j! ) +w aw) 

where 


a = 



Vi = 7T 


(fa - fa) 


log ( 


|x o + (a,0)| 
|x o -(a,0)| 


V 2 = 0.5 (fa + fa) 


(117) 
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with xo = (|). For each value of d, we compute the charge q\ (since q -2 = —qi), 
the iteration count for the elastance problem nn,eU and the relative L 2 error 


of the potential on boundary Fj given by e, = 


f r . l“-Mex | 2 


We emphasize 


f r . IU e x\ 2 d. 

again that this is not just a test of backward stability for the elastance solver, 
since we are solving two different boundary value problems. 


d 

Qi 

flit,el 

ei 

e-2 

0.5 

-0.239487 

4 

5.910 -8 

1.5 10 _Y 

0.05 

-0.743917 

8 

2.010" 5 

3.3 10 _t> 

0.005 

-2.348079 

15 

3.310 -b 

5.110 _t> 


Table 1: Summary of results for the capacitance and elastance problems with 
two discs. 

4.1.2 Splash test 

We repeat the test above with a more complicated geometry. We now consider 
5 conductors Dj , whose boundaries F ; - are parametrized by 

Xj (9) = Xj + rj(9) cos (9 + /3j) (118) 

Vj (°) = Vj + r j(0) sin ( 0 + Pj) (H9) 

where 

12 

rm = i+E a,j t k sin (k6), (120) 

fc =i 

with the coefficients cij t k are uniformly chosen from [0,0.1] and prescribe an 
arbitrary potential on each of these objects. 

We list here the parameters for defining the geometry and the exact solution 
in the previous section. The table of centers ir'j, j/j : , and fd 3 is given below. 



Ti 

r 2 

r 3 

r 4 

r 5 


- 1.2 

1.2 

0 

- 1.2 

1.2 


0 

0 

- 2.2 

-4.4 

-4.4 

Pi 

7 T 

0 

7T 

_8_ 

37T 

4 

7r 

4 


Table 2: Parameters for setting up splash test for the Elastance and Mobility 
problems. 

In the next table, we list the coefficients a 3t k for j = 1,2... 5 and k = 

1 , 2 ,... 12 . 
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Fi 

r 2 

r 3 

r 4 

r 5 

0.012065 

0.017038 

0.070082 

0.029959 

0.012613 

0.064385 

0.041668 

0.094629 

0.069290 

0.004017 

0.006234 

0.011991 

0.046520 

0.005102 

0.07413 

0.049028 

0.022743 

0.038905 

0.067634 

0.052361 

0.030608 

0.035266 

0.043884 

0.089215 

0.084973 

0.081641 

0.10864 

0.030143 

0.097489 

0.002916 

0.099718 

0.087338 

0.084480 

0.004693 

0.081962 

0.042460 

0.096291 

0.008018 

0.055024 

0.020443 

0.076748 

0.053323 

0.069852 

0.085238 

0.069016 

0.084684 

0.040564 

0.047617 

0.070539 

0.056950 

0.016811 

0.085034 

0.015078 

0.069771 

0.051020 

0.040454 

0.016044 

0.050553 

0.051137 

0.092286 


Table 3: Coefficients (ij,k- For fixed j, the coefficients a : j,k for T,- are listed in 
order of increasing k. 


The prescribed potentials (f>j on T ? are given below, as well as the L 2 norms 

/ Jr |v — U e x | 2 ds x 

for the errors in it, given by = * —j —,- ^ - on the boundary T^. The 

potential u is computed after solving the elastance problem to see if we recover 
the exact values u ex |r = 4>i- 


j 

1 

2 

3 

4 

5 

(f>j 

0.120625 

0.643859 

0.062342 

0.490279 

0.306079 

e o 

2.110" 5 

4.2 10” B 

2.4 10 -b 

8.2 10~ B 

8.010 -e 


Table 4: Prescribed potential on the boundary and the relative L 2 error in 
potential on the boundary. 

The elastance problem converged in 30 GMRES iterations, with a relative 
residual of 10 -6 . In Fig. [5j we show a contour plot of it, with boundary values 
set to the 4>j. 
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Figure 5: Contour plot of the potential u in the exterior of U jDj. 


4.1.3 Application: Computing dielectric properties of nanocompos¬ 
ites 


Nanoconrposites are composite material consisting of nanoparticles in a host 
medium. Of particular interest are nanocomposites consisting of metallic parti¬ 
cles in a homogeneous organic host due to their applications in transformation 
optics and high energy density storage materials. We shall treat the nanocom¬ 
posite as a collection of nanoparticles which are perfect conductors in ambient 
space. Computing bulk dielectric properties of such materials as a function of 
shape, orientation and the volume fraction of these nanoparticles is of practical 
interest. Low frequency dielectric constants are typically determined experimen¬ 
tally using “capacitance” measurements. The dielectric constant is determined 
by measuring the voltage drop between two charged plates in the presence and 
absence of the nanocomposite. If the two conducting plates have charge ±Q 
and the measured potential difference is AV, then the “capacitance” of the 
configuration is computed as 



( 121 ) 


The potential drop, AV can be computed by solving an elastance problem. 


Remark 8. It should be noted that obtaingin C in this manner is different 
from computing the mutual capacitance between the two plates for the given 
configuration of nanoparticles fffff] /. Experimentally one could have applied a 
potential difference between the two plates and measured the charge accumulated 
on them. However, to determine the mutual capacitance of this configuration 
numerically, one would need to know the potentials on each of the nanoparticles, 
and this data is not available. 


For fixed volume fraction 
study in 


we carry out a two-dimensional version of the 
In particular, we study the effects of varying the number of 
particles and their aspect ratio. Let D\ and D 2 with boundaries Ti and T 2 , 
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represent the capacitor plates. The boundaries Ti and T 2 are shifted copies of 
a rounded bar 7 parametrized by 



ioos 2 _|_ g . erf (10s)^ 

»6( f,f] 




( 122 ) 

(123) 


The curve 7 is discretized by sampling it at Sk = — § + n(k — 0.5 )/N, k = 
1,2,... 2 N. We verify that the curve is well-resolved by studying the discrete 
Fourier coefficients of the sampled curve and choose N sufficiently large that 
the curve is approximated to at least the desired tolerance. More precisely, the 
boundaries Ti and T 2 are parametrized by {x (s), y (s) ± 1 . 1 ) and discretized 
using 800 points each. For the nanoparticles, we use an m x 10 lattice of elliptic 
inclusions, each of which has an area equal to and an aspect ratio A. They 
are centered at 


-0.9 + 


1.8 (k — 1 ) 

1CT~ 


,-0.9 + 


1-8 (j — 1 ) 


j = 1,2. . .m, k = 1,2... 10. (124) 


The aspect ratio A is restricted to ensure that the nanoparticles do not overlap 
Each of these elliptical inclusions is discretized using an equispaced sampling of 
the central angle with 600 points each. 


C 




C 




Figure 6 : Capacitor plates with intervening nanoparticles. Here, there are m = 4 
rows with aspect ratio set to A = 0.5. 


In this case, we have a total of 10m + 2 conductors whose boundaries are 
discretzed with N pts = 6000m + 1600 points. We prescribe charges 1 and —1 on 
conductors Ti and T 2 , respectively, and assume the other 10 m conductors are 
charge neutral. We measure the potential difference AV = V\ — V 2 where V) is 
the potential on F^, i = 1,2 and compute the capacitance via equation (1211 
for various values of m, n and A. As before, we compute the layer potentials 
using a 6 th order GLQBX scheme accelerated via an FMM, and iterate using 
GMRES until the relative residual in our computation is less than 10 -6 . 

Let Cq be the capacitance in the absense of the nanocomposite (corresponding to 
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m 

A 

Npts 

N it 

^solve 

C 

0 

- 

1600 

8 

0.4539 

2.2949 

1 

0.25 

7600 

11 

3.2965 

2.3147 

0.5 

8 

2.2417 

2.3073 

1.0 

8 

2.2907 

2.3033 

2.0 

8 

2.2167 

2.3013 

4.0 

11 

3.4305 

2.3003 

4 

0.25 

25600 

11 

11.5612 

2.3191 

0.5 

9 

8.7257 

2.3095 

1.0 

7 

7.2139 

2.3047 

2.0 

9 

8.6497 

2.3023 

4.0 

10 

10.1145 

2.3012 


0.25 


11 

44.8332 

2.3200 


0.5 


9 

33.6639 

2.3099 

16 

1.0 

97600 

8 

31.2122 

2.3049 


2.0 


9 

33.5139 

2.3024 


4.0 


11 

46.4239 

2.3012 


Table 5: Capacitance of nanocomposites. Nu is the number of GMRES itera¬ 
tions, and t so i ve is the time taken to solve the elastance problem in seconds on 
a single CPU core. 


to = 0 ). We plot below the percentage change in capacitance 100 (c — Co) /Co 
as a function of m and A. 


Aspect Ratio (A) 


m=l 

m=4 

m=16 


Figure 7: % change in C as a function of to and A. 


4.2 Mobility problem 

We turn now to a test for our mobility representation. Given prescribed veloc¬ 
ities for a set of rigid bodies, we solve the resistance problem and compute the 
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resulting forces and torques on them. We then use these forces and torques as 
input for the mobility problem and check that the velocity on the boundary of 
the rigid body is the prescribed rigid body motion. As before, this is a stringent 
test, since the integral equation for the resistance problem is not simply the 
inverse of the integral equation for the mobility problem. 

We consider, as above, the domain exterior to N rigid bodies whose 
boundaries are given by IV We prescribe velocities u = v* + w* (x — x^) 1 " on 
and, solve the resistance problem to compute the forces and torques on the 
rigid bodies Di and also the velocity at oo, Uoo = limi*^,*, u (x). We use 
these forces and torques as input for the mobility problem to compute the rigid 
body motions. Let p t mob denote the incident velocity field due to the forces 
and torques as described in section 3.2 and let p i mob represent the unknown 
density on r, for the mobility problem. To summarize, we set 


u (x) = u inc (x) + u sc (x) + Uoo = S r (p mob (x) + p mob (x)) + , (125) 

and wish to solve 


-I + K, + L ) p mob — - ( -I + K. ) p mob 


(126) 


where p mob = (n 1>mob , p 2 ,mob) - Pmob = (Pi,mob, P2,mob) , u oo is the velocity at 
oo computed in the resistance problem, and the operators K. and L are described 
in section |3.5[ This is a small modification of the representation presented in 
to account for the velocity at oo. After solving for p mob , we verify 


section 3.2 


that u in equation (125) is the original rigid body motion. 


4.2.1 Two discs 

We hrst test the mobility representation in the exterior of two discs, with the 
same geometry is above, in section [TTTT] However, we use a finer discretization 
to resolve the more singular densities incurred in the mobility problem. 



Figure 8: Discretization of the discs for Mobility example. 
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Refering to Fig. [sj we set u|r 4 = v* + Wi (x — for i = 1,2 (i = 
corresponds to the disc on the left), where we set Vi = (2.09,1.00), v 2 = 
(—1.034,0.254), to i = 0.12 and w 2 = 0.33. We again test the problem for 


d = 0.5,0.05,0.005. 

As above, we use a Nystrom discretization with Gauss-Legendre panels, 
s i t j t i denotes the j th Gauss-Legendre node on panel i on boundary l. Let 


Pi,j,l,mob^i,j,l,mob denote the densities at s 
ture scheme for evaluating the layer potential, 


We use the GLQBX quadra- 
I + K, in equation (126). We 


use an iterative GMRES-based solver to obtain p mob , with a relative residual 
tolerance of 10~ 6 . The physical conditioning increases as d —> 0, requiring a 
large number of iterations. To impr ove the rate of convergence of GMRES, we 
use L 2 weighting for the unknowns 
as unknowns. 
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. That is, we use = Hi )jt i trnoby /n 


We solve the following linear system 


Dl-I + JC 


L D 




K Pn 


',ob 5 


(127) 


where K. and L are discretized versions of /C and L respectively, and D is the 
diagonal operator given by Dn mob = 

We plot below the net surface traction p 1 mob + Pi mob and a quiver plot of 
the velocity field in the exterior of the two discs for d = 0.05 (Figs. [9] and 10). 
From symmetry considerations, p 2 , mob + p 2 ,mob = - (Pi,mob + Pi,mlb)- 




Figure 9: Integral equation solution for the mobility problem: pi t i, mo b+Pi,i,mob 
(left) and p 2 ,i, mo b + p 2 ,i, mob (right) as a function of 9 for d = 0.5, 0.05, 0.005 
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2.4 


1.6 



Figure 10: Quiver plot of u and contour plot of |u| in the exterior of the two 
discs. 


For each value of d , we compute the forces and torques Fi,i, T\, T 2 
(since F 2 = —FQ, the iteration count for the mobility problem, tin, and the 
relative L 2 error of the velocity on both the boundaries Ti and T2 given by 

I f r |u— u ex \ 2 ds x 

-4 —:-T 2-7 -. This is again a nontrivial test of our solvers, as we 

J p. I Hex| 

enforce boundary conditions on the fluid stress here, not the velocity u. 


d 

Fi, 1 

t*2,l 

Ti 

t 2 

nu 

ei 

e2 

0.5 

27.180434 

-6.575686 

-1.496082 

1.494675 

7 

8.810-* 

1 . 810- 5 

0.05 

499.08688 

-15.202716 

-11.159661 

-4.859692 

19 

5.110~ B 

8.5 lO - * 

0.005 

14653.544 

-40.877338 

-42.867299 

-24.078713 

60 

1.010" B 

2.2 10“ b 


Table 6 : Summary of results for two discs test for resistance and mobility prob¬ 
lems. 


4.2.2 Splash test 

We repeat the test above in a more complicated geometry, but consider 5 bodies 
Dj , whose boundaries Tj are parametrized as 


Xj ( 6 ) = Xj + vj{0) cos (9 + ftj) 

Vj ( 0 ) = V c j + r j(9) sin ( 0 + Z 3 .?) 


where 


12 

rj(9) = l + E cij,k sin ( k0 ) 
k =1 


(128) 

(129) 


( 130 ) 


where the parameters ay*,, x? and y c , are the described in section 


4.1.2 
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The prescribed velocities vj = (vij,V 2 j) and cuj on Tj are given below, 

o . . . / fr l u—Uea; 1 2 ^ s x 

along with the Lr norm in the error in u, given by = w — f —;--on the 

V l Uea: I dSx 

boundary Tj, after solving the mobility problem is listed below where Uexliq = 
Vi + Uli (x - X?) . 


j 

1 

2 

3 

4 

5 

v l,i 

-0.379375 

-0.009720 

0.497180 

0.346837 

-0.197527 


0.143846 

-0.193921 

-0.075401 

-0.331891 

0.273004 

COj 

-0.437658 

0.316414 

0.267477 

-0.095456 

-0.184353 

e 7 

1.810" 5 

2.510” b 

1 .110“ 5 

1.310~ 5 

1.310~ 5 


Table 7: Prescribed velocity on the boundary and the relative L 2 error in velocity 
on the boundary after solving the resistane and mobility problems. 


The mobilitiy problem converged in 71 GMRES iterations to a relative tol¬ 
erance of 10” 6 . We show a quiver plot for the velocity u corresponding to the 
above prescribed value of velocity on the boundary u|r - ■ The background is a 
contour plot of the magnitude |u| in the exterior of the U jDj (Fig. 111. 




1.2 


0.8 



Figure 11: Quiver plot of u superimposed on contour plot of |u| in the exterior 
of U jDj. 


5 Conclusions 


We have derived a new, physically motivated integral formulation for the elas- 
tance problem in exterior domains. The analogous physical reasoning yields a 
new derivation of an integral equation developed earlier by Kim and Karrila 
for the mobility problem 27 . Discretization of the resulting integral equations 
using the quadrature scheme GLQBX 33 34 permits high order accuracy to 
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be obtained in complex geometry, including the interaction of close-to-touching 
boundary components. The resulting linear systems can be solved iteratively 
using GMRES and the necessary matrix-vector multiplications can be acceler¬ 
ated using the fast multipole method. If N denotes the number of points used 
in the discretization of the physical boundaries, the total cost scales linearly 
with N. We are currently working on the extension of our scheme to closed or 
periodic systems and to problems in three dimensions. 


Acknowledgments 

This work was supported in part by the Applied Mathematics program in the 
Department of Energy, Office of Advanced Scientific Computing Research, under 
contract DEFGO288ER25053 and by the Office of the Assistant Secretary of 
Defense for Research and Engineering and AFOSR under NSSEFF Program 
Award FA9550-10-1-0180. The authors would like to thank Alex Barnett and 
Mike O’Neil for several useful discussions. 


References 

[1] W. R. Smythe, Static and Dynamic Electricity. McGraw Hill, New York, 
1975. 

[2] J. D. Jackson, Classical Electrodynamics. Wiley, New York, 1975. 

[3] S. Kapur and D. E. Long, “Large-scale capacitance calculation,” in in Proc. 
37th Design Automation Conference, pp. 744-749, 2000. 

[4] Y. C. Pan, L. X. Wan, and W. C. Chew, “A fast multipole method based 
calculation of the capacitance matrix in a stratified medium,” in Anten¬ 
nas and Propagation Society International Symposium, 2000. IEEE, vol. 4, 
pp. 1876-1879 vol.4, July 2000. 

[5] J. Tausch and J. White, “Second-kind integral formulations of the capac¬ 
itance problem,” Advances in Computational Mathematics, vol. 9, no. 2, 
pp. 217-232, 1998. 

[6] W. T. Weeks, “Calculation of coefficients of capacitance of multiconduc¬ 
tor transmission lines in the presence of a dielectric interface,” Microwave 
Theory and Techniques, IEEE Transactions on, vol. 18, pp. 35-43, Jan 

1970. 

[7] A. M. Chang and J. C. Chen, “The Kondo effect in coupled-quantum dots,” 
Reports on Progress in Physics, vol. 72, p. 096501, Sept. 2009. 

[8] W. G. van der Wiel, S. De Franceschi, J. M. Elzerman, T. Fujisawa, 
S. Tarucha, and L. P. Kouwenhoven, “Electron transport through double 
quantum dots,” Reviews of Modern Physics, vol. 75, pp. 1-22, 2003. 


35 



[9] J. Koch, T. Yu, J. Gambetta, A. Houck, D. Schuster, J. Majer, A. Blais, 
M. Devoret, S. Girvin, and R. Schoelkopf, “Charge-insensitive qubit design 
derived from the cooper pair box,” Physical Review A, vol. 76, p. 042319, 
Oct. 2007. 

[10] S. Delong, F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, 
and A. Donev, “Brownian Dynamics without Green’s Functions,” 
arXiv: lfOl-4198v2, 2014. 

[11] K. Icliiki and J. F. Brady, “Many-body effects and matrix inversion in low- 
Reynolds-number hydrodynamics,” Physics of Fluids, vol. 13, no. 1, p. 350, 
2001 . 

[12] S. Kim and S. Karrila, Microhydrodynamics: Principles and Selected Ap¬ 
plications. Butterworth - Heinemann series in chemical engineering, Dover 
Publications, 2005. 

[13] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized 
Viscous Flow. Cambridge Texts in Applied Mathematics, Cambridge Uni¬ 
versity Press, 1992. 

[14] W. C. Chew, J.-M. Jin, E. Michielssen, and J. Song, Fast and Efficient Al¬ 
gorithms in Computational Electromagnetics. Artech House, Boston, 2001. 

[15] Y. Liu, Fast Multipole Boundary Element Method: Theory and Applications 
in Engineering. Cambridge University Press, New York, 2009. 

[16] L. Greengard and V. Rokhlin, “A new version of the fast multipole method 
for the Laplace equation in three-dimensions,” Acta Numerica , vol. 6, 
pp. 229-269, 1997. 

[17] N. Nishimura, “Fast multipole accelerated boundary integral equation 
methods,” Appl. Mech. Rev., vol. 55, pp. 299-324, 2002. 

[18] A. Greenbaum, L. Greengard, and G. McFadden, “Laplace’s equation and 
the Dirichlet-Neumann map in multiply connected domains,” Journal of 
Computational Physics, vol. 105, pp. 267-278, 1993. 

[19] J. Helsing and E. Wadbro, “Laplace’s equation and the Dirichlet-Neumann 
map: a new mode for Mikhlin’s method,” Journal of Computational 
Physics, vol. 202, pp. 391-410, Jan. 2005. 

[20] S. Mikhlin, Integral equations and their applications to certain problems: 
in mechanics, mathematical physics and technology. International series of 
monographs in pure and applied mathematics, Macmillan, 1964. 

[21] K. Nabors, F. T. Korsmeyer, F. T. Leighton, and J. White, “Pre¬ 
conditioned, adaptive, multipole-accelerated iterative methods for three- 
dimensional first-kind integral equations of potential theory,” SIAM Jour¬ 
nal on Scientific Computing, vol. 15, pp. 713-735, 1994. 


36 



[22] G. Biros, L. Ying, and D. Zorin, “A fast solver for the Stokes equations 
with distributed forces in complex geometries,” Journal of Computational 
Physics, vol. 193, pp. 317-348, 2003. 

[23] L. Greengard, M. C. Kropinski, and A. Mayo, “Integral equation methods 
for Stokes flow and isotropic elasticity in the plane,” Journal of Computa¬ 
tional Physics, vol. 125, pp. 403-414, 1993. 

[24] H. Power and G. Miranda, “Second kind integral equation formulation of 
Stokes flows past a particle of arbitrary shape,” SIAM Journal on Applied 
Mathematics, vol. 47, pp. 689-698, 1987. 

[25] R. Cortez, L. Fauci, and A. Medovikov, “The method of regularized 
Stokeslets in three dimensions: Analysis, validation, and application to 
helical swimming,” Physics of Fluids, vol. 17, no. 3, p. 031504, 2005. 

[26] M. Kropinski, “Integral equation methods for particle simulations in creep¬ 
ing flows,” Computers & Mathematics with Applications, vol. 38, pp. 67-87, 
Sept. 1999. 

[27] S. J. Karrila and S. Kim, “Integral equations of the second kind for Stokes 
flow: direct solution for physical variables and removal of inherent accuracy 
limitations,” Chemical engineering communications, vol. 82, no. 1, pp. 123- 
161, 1989. 

[28] L. Af Klinteberg and A.-K. Tornberg, “Fast Ewald summation for Stoke¬ 
sian particle suspensions,” International Journal for Numerical Methods in 
Fluids, vol. 76, no. 10, pp. 669-698, 2014. 

[29] R. Kress, Linear Integral Equations. Applied Mathematical Sciences, 
Springer New York, 1999. 

[30] R. Guenther and J. Lee, Partial Differential Equations of Mathematical 
Physics and Integral Equations. Prentice Hall, 1988. 

[31] J. Sifuentes, Z. Gimbutas, and L. Greengard, “Randomized methods for 
rank-deficient linear systems,” arXiv:lf01.3068, 2014. 

[32] R. Finn and W. Noll, “On the uniqueness and non-existence of Stokes 
flows,” Archive for Rational Mechanics and Analysis, vol. 1, no. 1, pp. 97 
106, 1957. 

[33] M. Rachh, Integral equation methods for problems in electrostatics, elasto- 
statics and viscous flow. 2015. Thesis (Pli.D.)-New York University. 

[34] M. Rachh, A. Klockner, and M. O’Neil, “A hybrid version of QBX (quadra¬ 
ture by expansion) for the evaluation of layer potentials in close-to-touching 
geometries,” in preparation. 


37 



[35] A. Klockner, A. Barnett, L. Greengard, and M. O’Neil, “Quadrature by 
expansion: A new method for the evaluation of layer potentials,” Journal 
of Computational Physics , vol. 252, pp. 332-349, Nov. 2013. 

[36] J. Bremer and V. Rokhlin, “Efficient discretization of Laplace boundary in¬ 
tegral equations on polygonal domains,” Journal of Computational Physics, 
vol. 229, no. 7, pp. 2507 2525, 2010. 

[37] A. Barnett, “Evaluation of layer potentials close to the boundary for 
Laplace and Helmholtz problems on analytic planar domains,” SIAM J. 
Sci. Comput., vol. 36, pp. A427-A451, 2014. 

[38] X. Zheng, J. Fontana, M. Pevnyi, M. Ignatenko, S. Wang, R. Vaia, and 
P. Palffy-Muhoray, “The effects of nanoparticle shape and orientation on 
the low frequency dielectric properties of nanocomposites,” Journal of Ma¬ 
terials Science , vol. 47, no. 12, pp. 4914-4920, 2012. 


38 



